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Abstract 

We use equivariant bifurcation theory to investigate pattern selection at the 
onset of a Turing instability in a general two-component reaction-diffusion sys- 
tem. The analysis is restricted to patterns that periodically tile the plane in 
either a square or hexagonal fashion. Both simple periodic patterns (stripes, 
squares, hexagons, and rhombs) and "superlattice" patterns are considered. The 
latter correspond to patterns that have structure on two disparate length scales; 
the short length scale is dictated by the critical wavenumber from linear theory, 
while the periodicity of the pattern is on a larger scale. Analytic expressions 
for the coefficients of the leading nonlinear terms in the bifurcation equations 
are computed from the general reaction-diffusion system using perturbation the- 
ory. We show that, no matter how complicated the reaction kinetics might be, 
the nonlinear reaction terms enter the analysis through just four parameters. 
Moreover, for hexagonal problems, all patterns bifurcate unstably unless a par- 
ticular degeneracy condition is satisfied, and at this degeneracy we find that the 
number of effective system parameters drops to two, allowing a complete char- 
acterization of the possible bifurcation results at this degeneracy. For example, 
we find that rhombs, squares and superlattice patterns always bifurcate unstably. 
We apply these general results to some specific model equations, including the 
Lengyel- Epstein CIMA model, to investigate the relative stability of patterns as a 
function of system parameters, and to numerically test the analytical predictions. 



1 Introduction 

Regular patterns arise in a wide variety of physical, chemical and biological systems that 
are driven from equilibrium. The origin of these spatial patterns can often be traced to 
a symmetry-breaking instability of a spatially-homogeneous state. Such an instability 
mechanism was proposed by Alan Turing in 1952 |l| for pattern formation in chemical 
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systems featuring both reaction and diffusion. The key insight of Turing was that 
diffusion, which is usually thought of as a stabilizing mechanism, can actually act as a 
destabilizing mechanism of the spatially-uniform state. This is the case, for example, 
in a class of two-component activator-inhibitor systems, provided the inhibitor diffuses 
more rapidly than the activator. This competition between long-range inhibition and 
short-range activation can lead to spatially inhomogeneous steady states called Turing 
patterns j^]. Experimental evidence of Turing patterns was discovered only recently]^, 
and a number of experiments on the Turing instability have now been carried out on 
the chlorite- iodide-malonic acid (CIMA) reaction in a gel § ||. These experiments 
have demonstrated the existence and stability of a variety of regular patterns, including 
stripes, hexagons, and rhombs. At the same time, a two variable model of the CIMA 
reaction has been proposed by Lengyel and Epstein . 

Equivariant bifurcation theory is a powerful tool for analyzing the nonlinear evo- 
lution of symmetry-breaking instabilities in pattern-forming systems. This approach 
determines the generic behaviors associated with bifurcation problems that are equiv- 
ariant with respect to a given group of linear transformations. In particular, the form 
of the bifurcation problem is restricted by symmetry considerations, limiting the be- 
havior to some finite set of possibilities. The problem of determining which type of 
behavior occurs in a given physical system then reduces to one of computing the co- 
efficients of specific nonlinear terms in the bifurcation problem. In this paper we 
perform this computation for a general two-component reaction-diffusion system in a 
neighborhood of a Turing bifurcation. We consider two different symmetry groups for 
the problem, one associated with Turing patterns that tile the plane in a square lat- 
tice, and the other associated with patterns that tile the plane in a hexagonal lattice. 
The symmetry groups are D4+T^ and Dg-i-T^, respectively, where D„ {n = 4,6) is 
the dihedral group of symmetries of the fundamental tile, and T^ is the two-torus of 
translation symmetries associated with the doubly-periodic solutions. The fundamen- 
tal irreducible representations of these groups, which are four- and six-dimensional, 
lead respectively to bifurcation problems that address competition between stripes and 
simple squares, and between stripes and simple hexagons. The higher dimensional ir- 
reducible representations, which are eight- and twelve-dimensional, lead to bifurcation 
problems that address competition between these simple Turing patterns and a class of 
patterns that are spatially-periodic on a larger scale than the simple patterns, but still 
arise in a primary bifurcation from the homogeneous state at the Turing bifurcation 
point H, ^, . An example of such a "superlattice" pattern was recently observed 
in an experiment on parametrically excited surface waves [p^ . Finally, the higher- 
dimensional group representations also allow us to investigate the relative stability 
properties of rhombs and squares or rhombs and hexagons. 

A number of earlier studies have used bifurcation theory to investigate Turing pat- 
tern formation for reaction-diffusion systems. For instance, Othmer described a 
group theoretic approach for analyzing the Turing bifurcation nearly twenty years ago, 
in 1979. Recently Callahan and Knobloch derived and analyzed the general equiv- 
ariant bifurcation problems for three-dimensional Turing patterns having the spatial 
periodicity of the simple cubic, face-centered cubic and body-centered cubic lattices. 
They applied their general results to the Brusselator and the Lengyel-Epstein CIMA 
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reaction-diffusion models Gunaratne, et al. developed an approach to pat- 

tern formation problems based on a model equation akin to Newell- Whitehead-Segel, 
but with full rotational symmetry. They used this model to investigate bistability of 
hexagons and rhombs, and tested their predictions against laboratory experiments. 

In addition to these general bifurcation studies of Turing pattern formation, a 
number of specific reaction-diffusion models have been investigated analytically and 
numerically. Dufiet and Boissonade investigated numerically the formation of 
two-dimensional Turing patterns for the Schnackenberg reaction-diffusion system. The 
formation of three-dimensional Turing structures has been investigated for the Brus- 
selator by DeWit, et al. jl^. Turing patterns for the two- variable reaction-kinetics 
based model of the CIMA reaction, proposed by Lengyel and Epstein have been 
investigated using bifurcation theory [|l9| , as well as numerically . 

This paper extends existing bifurcation studies of Turing patterns in reaction- 
diffusion systems in two primary ways. The first, as discussed above, is that we consider 
the higher dimensional (8 and 12, respectively) irreducible group representations asso- 
ciated with the square and hexagonal lattice tilings. This corresponds to choosing a 
fundamental tile for the periodic patterns large enough to admit many copies of the 
simple Turing patterns - stripes, hexagons, squares, and rhombs - as well as new solu- 
tions in the form of "superlattice" patterns. In this way we extend the standard stripes 
vs. hexagons stability analysis, while also determining conditions under which other 
more exotic Turing patterns might arise. 

The second way in which our analysis differs from earlier studies is that it does 
not restrict attention to a specific model {e.g. Schnackenberg, Brusselator, Lengyel- 
Epstein CIMA model, etc.). The coefficients in the bifurcation equations are derived for 
a general two-component reaction diffusion system. This allows us to make a number 
of general statements about Turing patterns in these systems. For example, we show 
that the bifurcation results depend on four effective parameters that may be simply 
computed from the nonlinear terms in the reaction-kinetics. Moreover, in the case of 
the standard degenerate hexagonal bifurcation problem, we show that the details of the 
reaction-kinetics enter the analysis through just two parameters. Even more surprising 
is the result that at this degenerate point all angle-dependence, as associated with 
rhombic patterns, drops out of the corresponding bifurcation problems; in this case 
we show that all rhombic patterns are unstable to rolls, even when hexagons are the 
preferred planform. This is in contrast to the general stability analysis of Gunaratne 
et al. ||l6|, which suggest that rhombs with angle close to 7r/3 should be stable when 
hexagons are stable. Thus our analysis indicates that two-component reaction-diffusion 
models may have a more rigid mathematical structure than anticipated. 

Our paper is organized as follows: Section ^ reviews the necessary conditions for 
a Turing instability in a two-component reaction-diffusion system. It also provides 
the necessary background on the mathematical framework for the bifurcation analysis. 
Section ^ calculates the bifurcation equation coefficients using perturbation methods, 
with the final results summarized in an Appendix. Section ^ describes the concept of 
parameter collapse, where multiple system parameters collapse into a small number 
of effective parameters. Special attention is given to the degenerate bifurcation prob- 
lem on the hexagonal lattice, where we show that the details of the reaction-kinetics 
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enter the analysis in this case through only two effective parameters. The resulting 
predictions are then then tested numerically; specifically, we present an example of 
tri-stability between stripes, simple hexagons, and super hexagons in a neighborhood 
of the degenerate bifurcation problem. Section || applies the general results to a model 
of the CIMA reaction proposed by Lengyel and Epstein Q and a model demonstrating 
super squares. Finally, we suggest some directions for future work and summarize the 
key points of the paper in the conclusions section. 



2 Formulation of Bifurcation Problem 

2.1 Linearized Problem 

We begin by summarizing the conditions for a Turing instability to occur as some 
parameter is varied in the reaction-diffusion system 

Ut = V^u + f{u,v) 

Vt = K\/^v + g{u,v), "V^^d^^ + dyy. (1) 

We assume that the species diffuse at different rates and that v measures the concen- 
tration of the more rapidly diffusing species. Thus the ratio of diffusion coefficients 
is K > 1. The Turing instability first occurs as a symmetry-breaking steady state 
bifurcation of a spatially-uniform steady state of (|l|) ; associated with the instability is 
a critical wave number Qc 0, which is determined below. 

Without loss of generality, we take the spatially-uniform state to be u = u = 0. We 
expand f{u,v) and g{u,v) about this state as 

f{u,v) — au + bv + . . . 
g{u,v) — cu + dv + . . . 

yielding the linearized reaction-diffusion system 

ut\ fa + \7^ b \ f u\ 



Vt J \ c d + KV^ J \ v ^ 
The stability of the spatially uniform solution is determined by substituting 



into (H) to obtain the eigenvalue problem 



<j{q)^ = L{q)(; L{q)^(^^ f d-Kq')' 

The eigenvalues (Ti(g), (J2{q) depend on the parameters a,b, c, d, K , as well as the 
wavenumber of the perturbation q. A Turing bifurcation occurs for parameter values 
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where there is a zero eigenvalue for some q = qc ^ Q with all other modes being damped, 
i.e., 

(Tl((?c) < 0-2 (gc) = 0, 
Re(crj (g)) < 0, j = 1, 2, for all q ^ q^. 
These inequalities imply that 

Tr(L(<j)) < 0, Vq, (4) 
Det(L(g)) > 0, q^qc, 

at the bifurcation point. The critical wave number q^ is given by 

Det(L(<z,)) = 0. (5) 

Since Det{L{q)) is quadratic in q^ and concave up, the second of conditions (^) is 
assured if qc is a minimum of Det{L{q)): 



dT)et{L{q)) 



dq 



= 0. (6) 

9=?c 



Thus the parameters at the Turing bifurcation (a = ttc, b = be, c = Cc, d = dc, K 
Kc, q = qc) satisfy 



2 Kcac + dc 



(Kcttc - dcf + AKcbcCc = 0, (8) 

and must also satisfy the inequalities 

KcUc + dc > 0, 

bcCc < 0, (9) 
ac + dc < 0, 

Qcdc - bcCc > 0, 

to ensure that q^ > and that perturbations at g = decay. From the initial as- 
sumption that K > 1, it follows that Gc > and dc < 0, which for be < 0, Cc > is 
consistent with the interpretation of u as the concentration of the activator and v as 
the concentration of the inhibitor in (|l|). 

2.2 Spatially-Periodic Solutions of the Nonlinear Problem 

For the two-dimensional spatially-unbounded problem there are infinitely many neu- 
trally stable Fourier modes at the bifurcation point: all Fourier modes associated with 
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wave vectors on the critical circle |q| = qc are neutrally stable. In the remainder of the 
paper we focus on spatially doubly-periodic solutions of (|l|): 

u{y, t) = u{t + di, i) = u{y + d2, t), i;(r, t) = v{y + di, t) = u(r + d2, t), 

where r = (x, y) and di, d2 are linearly independent vectors in R^. We further restrict 
to the cases where |di| = |d2|. We then write each field u,v in terms of its Fourier 
series, e.g. 

u{v,t) = ^ ^i,,fe(t)e'»(^"'^^+'='^^)- , u.,..k^ul^, (10) 

where ki and k2 are unit vectors that satisfy 

gki • dj = 2-K5ij 

for an appropriate choice of q. Since the Fourier transform of spatially-periodic solu- 
tions is discrete, only a finite number of Fourier modes can have wave vectors q that lie 
on the critical circle |q| — qc at the bifurcation point; all other modes will be linearly 
damped. Thus, in this setting, finite-dimensional bifurcation theory may be applied in 
a neighborhood of the Turing bifurcation point. 

For Turing patterns that tile the plane in a square lattice the unit vectors ki and k2 
are separated by an angle of 7r/2, whereas in the hexagonal case they are separated by 
27r/3; all angles that are not a multiple of 7r/3 or 7r/2 correspond to rhombic lattices. 
The parameter q in (0), which is determined by the imposed period, also dictates the 
the spacing of the lattice points in the Fourier transform, and may in turn be used to 
determine which lattice points (if any) lie on the critical circle. Of particular interest 
here is the ratio of this scale to the natural scale of the problem q^. For instance, if 
q/qc ~ 1 then on the square lattice there arc four lattice points which lie on the critical 
circle: (j, fc) = (±1,0), (0, ±1) in ([l^), corresponding to Fourier modes 

If, however, a tighter spacing with q/qc = l/'s/S is used then there are eight square- 
lattice modes which lie on the circle: (j, fc) = (2,±1), (— 2,±1), (1,±2) and (— 1,±2). 
Examples for these cases, as well as analogous ones for the hexagonal lattice, are 
depicted in Figure |l|. 

In this paper we consider a discrete set of values of q, parameterized by an integer 
pair {m,n), which admit wave vectors in ( p^ with length qc, i.e., for which |g(mki + 
nk2)| = qc- Specifically, on the square lattice the values of q are related to {m,n) by 

q^=ql/{m^ + n^), (11) 

and on the hexagonal lattice by 

q^ = q^/ [m^ + — rrm) . (12) 
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Figure 1: Square and hexagonal lattices in Fourier-space plotted together with two 
examples of critical circles of radius qc. The dashed circles corresponds to Qc = q in 
( pi] ) and (12), i.e. {m,n) = (1,0); the solid circle corresponds to (171,71) ~ (2, 1) on the 
square lattice and (m, n) — (3, 1) on the hexagonal lattice. Note that the solid circle 
has twice as many critical modes as the dashed circle. 



In a neighborhood of the bifurcation point, {a,b, c, d, K) — (ac,bc,Cc,dc,Kc), we 
apply our analysis separately to the cases where there are eight critical Fourier modes 
on the square lattice and twelve on the hexagonal lattice (see Figure ||) . Specifically, 
on the square lattice, we investigate primary steady solution branches of the form 

zie'"^'-" + Z2e*«^^-"" + zse'"^'" + Z4e"?^^-"" + c.c. + • • • (13) 

where zi, . . . , Z4 are the complex- valued amplitudes of the critical modes, • • • represents 
the linearly-damped Fourier modes that are slaved to the critical modes, and Ki = 
(m, n), K2 = (— n,rn), K3 = (n, m), and K4 = (— m,n). Note that, while Ki • K2 = 
K3-K4 = 0, the angle between Ki and K3 depends on the integer pair (m, n). Similarly, 
on the hexagonal lattice. 



^ z^e'"^"-'' + C.C. + ■ ■ ■ (14) 



whereKi = m(l, 0)-Hn(-l/2, \/3/2) and K4 = m(l,0) + (m-n)(-l/2, V3/2). K2 and 
K3 are obtained by rotating Ki by ±27r/3, and K5 and Kg are obtained by rotating 
K4 by ±27r/3 (see Figure llj). Note that the angle between Ki and K4 again depends 
on (m, 71). 

Dionne and Golubitsky B used the equivariant branching lemma to prove that 
for the square lattice cases (|l^) , the following types of steady solutions arise in pitchfork 
bifurcations from the spatially-uniform state: 

• Stripes: z = (zi, Z2, z^, Z4) — (x, 0,0,0) 

• Simple squares: z — (x, x, 0, 0) 



7 



• Rhombs: z ~ (x, 0, x, 0) and (a;, 0, 0, x) 



• Super squares: z = {x,x,x,x) 

• Anti-squares: z — {x, x, —x, —x), 



where in each case x is a real amphtude. (Examples of super and anti-squares are given 
in Section 5^2, Figure ||). While the striped and simple square patterns are essentially 



the same for each m and n, the solutions in the form of rhombs, super squares and 
anti-squares differ for each co-prime m and n since the angle between Ki and K3 
(and Ki and K4) depends on m and n. The general form of the amplitude equations 
describing the branching and relative stability of the square lattice patterns is derived 
in 1^. To cubic order the bifurcation equations take the simple form 

Zl = iizi -f zi(ai|zip + a2\z2\'^ + aslzap -|- 04124^) H (15) 

with similar expressions for Z2, is, 24 and their conjugates. Here is the bifurcation 
parameter and ai , . . . , 04 are real. Linearizing (|l^) about a particular solution generates 
a set of eigenvalues which determines the stability of that solution with respect to 
perturbations on the lattice. The signs of these eigenvalues are given in Table m 
terms of ai, . . . , 04; a positive eigenvalue indicates instability. Note that higher order 
terms are necessary to determine the relative stability of the super square and anti- 
square patterns [^. In the next section we compute the coefficients ai, . . . , 04 needed 
to evaluate the signs of the eigenvalues in terms of the angle between Ki and K3 and 
in terms of parameters in the general reaction diffusion system (|^). 

Dionnc and Golubitsky Q also considered the hexagonal lattice problem and used 
group theoretic methods to determine that there are primary branches of the form 

• Stripes: z = (zi, Z2, 2:3, Z4, Z5, zg) = (a^, 0, 0, 0, 0, 0) 

• Simple hexagons: z = (x, x, x, 0, 0, 0) 

• Rhombs: z = (x, 0, 0, x, 0, 0), (x, 0, 0, 0, x, 0), and (x, 0, 0, 0, 0, x) 

• Super hexagons: z = (x, x, x, x, x, x) . 

Recently, Silber and Proctor ||Tl[| showed that there is an additional primary solution 
branch that has triangular symmetry: 

• Super triangles: z — (z, z, z, z, z, z) z G C, 

where, near the bifurcation point, the argument of the complex amplitude z is ^ IT- 
The stripes and simple hexagons are the same for each integer pair (TO,n), while the 
other solutions change with the values of m and n. The general form of the bifurcation 
equations describing the evolution of the critical modes zi, . . . , zq is derived in ^ . The 
cubic truncation is 

Zl = ^zi +7z*z* +zi(6i|zip +62(|22p + kap) 

+bi\zi\^ + b5\z5\^ + be\ze\^) + --- (16) 
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Table 1: Signs of eigenvalues for primary bifurcation branches on the square lattice; 
oi, . . . , Q4 are coefficients in the bifurcation equation ([15[). See ||9| for more details. 



Planform 


Signs of non-zero eigenvalues 


Stripes 
z (x, 0,0,0) 


sgn(ai), sgn(a2-ai), sgn(a3 - ai), sgn(a4 - ai) 


Simple Squares 
z = (x, X, 0, 0) 


sgn(ai+a2), sgn(ai - 02), sgn(a3 + 04 - ai - 02) 


Rhombs 
z = (x, 0, X, 0) 


sgn(ai + 03), sgn(ai - 03), sgn(a2 + 04 - ai - 03) 


Rhombs 

z = (x, 0, 0, x) 


sgn(ai + 04), sgn(ai - 04), sgn(a2 + 03 - ai - 04) 


Super Squares 

Z ('^; 


sgn(ai + a2 + 03 + 04), sgn(ai + a2 - as - 04) 
sgn(ai - a2 + 03 - 04), sgn(ai - a2 - as + 04) 
sgn(Co), where Co = 0(x2(™+"-i)) 


Ant i- Squares 


same as super squares, except Co ^ ~Co 
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with similar expressions for Z2,.-.,ze. While the stripes and rhombs solutions arise 
from the uniform state in a pitchfork bifurcation, all other solutions bifurcate trans- 
critically. The relative stability of the primary patterns are summarized from in 
Table |^. Note that unless the quadratic coefficient 7 is zero, all branches bifurcate 
unstably. Thus we will consider a degenerate bifurcation problem for which 7 w 0. 
Also note that, as in the case of super squares and anti-squares, the relative stability 
of the super hexagons and super triangles is determined by terms higher than cubic in 

® 0- 



3 Reduction to Bifurcation Equations 

Although symmetry determines the form of the bifurcation equations, the particular 
system determines o„ (and 6„) in equations (p^ (and ([l6|)). Each coefficient a„ is 
associated with a particular amplitude Zn in ( |lq ) , and each of these amplitudes is in 
turn associated with a Fourier wave vector K„. In this section we use perturbation 
theory to derive the coefficients a„, 7, and 6„ from the full system (|l|). 

The coefficients ai, . . . , 04 and 61 (= Oi), 64, 65, feg Txiay all be calculated by consider- 
ing a single bifurcation problem involving just two critical Fourier modes; the remaining 
coefficients 7 and &2 are calculated by considering a bifurcation to hexagons. To under- 
stand the need for only two modes in the computation of ai, . . . , 04, 64, 65, be consider, 
for example, the square lattice, where the K2, K3, and K4 vectors are all related to the 
Ki vector by an angle 9. Let 9j denote the angle between Ki and K^-; then 62 = tt/2 
and 64 =11/2 + 63 where ^3 is related to the integer pair (m, n) by 

2mn 

COS03 = ^— ^ (17) 

(see Figure |l|). The coefficients aj are a function of the angle 9j: 

a,=h{0„), J =2,3,4- 
Similarly, on the hexagonal lattice, 

bj^h{0,), J =4,5,6 

where 65 = 64 + 27r/3 and Oq = 64 + 4t:/3, with 0^ determined by the integers {m,n). 
From symmetry considerations we know that h{0) = h^n — 9). A calculation with 
two critical modes will determine h{9), and hence determine all of the angle-dependent 
coefficients in the bifurcation equations. 

To compute h{9) we seek solutions which are periodic on a rhombic (or square) 
lattice; specifically, solutions of the form 

u;,h = ^ie*i-'"-f Z9e*«-'- + c.c., (18) 
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Table 2: Branching equations and signs of eigenvalues for primary bifurcation branches 
on the hexagonal lattice; 7, 61, ... , be are coefficients in the bifurcation equation (16). 
See & ITcIl for more details. 



Planform and branching equation 


Signs of non-zero eigenvalues 


Stripes 
z = (x, 0,0, 0,0,0) 
= + bix^ + 0{x^) 


sgn(6i), 

sgn(7X + (&2 - bi)x'^), sgn(-7X + (62 - bi)x'^), 
sgn(64-5i), sgn(65-fei), sgn(66 - 61) 


Simple Hexagons 
z = {x, X, X, 0, 0, 0) 
= + 7x2 + ^ 262)2;^ + 0(x*) 


sgn(7X + 2(61 + 262)2;^), sgn(-7X + (&i - 62)2;^) 
sgn(-7x + (64 + b^ + be - bi - 262 )x2) 
sgn(— 7x + 0(x'^)) 


Rhombs 
z = (x,0,0,a;,0,0) 
= + (61 + 64)2;^ + 0(a;5) 


sgn(6i + 64), sgn(6i~64), sgn(Ci), sgn(C2), 
where (1+^2 = (-26i - 264 + 2^2 + ^5 + be)x'^, 
C1C2 = -l^x^ + {bi + 64 - 62 - 65) (&i +b4-b2- be)x^ 


Rhombs 
z = (x,0,0,0,a;,0) 
+ (61 + 65)2;^ + 0(a;-^) 


sgn(6i + 65), sgn(6i-65), sgn(Ci), sgn(C2), 
where (1 + ^2 = (-26i - 265 + 2^2 + 64 + fo6)2;^ 
C1C2 = -7^2;2 + (61 + 65 _ 62 - 64) (&i + fos - 62 - ^6)2;'' 


Rhombs 
z = (x, 0,0,0,0,2;) 
= Ai2; + (61 + 65)2;^ + 0(x5) 


sgn(6i + 66), sgn(6i-66), sgn(Ci), sgn(C2), 
where Ci + (2 = (-26i - 266 + 262 + 64 + &5)a;2 
C1C2 = -7^2;2 + (61 +be-b2- 64) (&i + ^6 - 62 - ^5)2;'' 


Super Hexagons 
= ^x + 7^2 + (&i + 252)x^ 

+ (64 + fo5 + 66).T^ + 0(x4) 


sgn(7x + 2(&i + 262 + 64 + &5 + be)x'^) 
sgn(7x + 2(61 + 2^2 - 64 - 65 - b6)x'^) 
sgn(-7x + 0(x3)), sgn(Ci), sgn(C2), 
where Ci + C2 = — 47X + 4(6i — 52)2;^ 
C1C2 - 4(7x - (61 - b2)x^)^ 

-2((64 - b,)^ + (64 - bef + (&5 - 66)'))2:4 
sgn(Co), where Co = 0(x2("-i)) 


Super Triangles 
z = (z,z,z,z,z,z), 
z = xe*'^,V'7^0,7r,... 

= ^iZ + JZ^ + {bl+2b2)\z\^Z 

+ (64 + 65 + fe6)|2|2z + 0(x4) 


Same as super hexagons 
except Co -Co 



11 



whereki — qc{l , 0) , kg = (cos 61, sin 0), and is not a multiple of 7r/3. The amplitudes 
evolve in a neighborhood of the bifurcation point /i = according to 

zi = ^j,zi + zi{ai\zi\'^ + h{d)\zef) 

ze = ^^zg + zgial\zg\^ + h{9)\zi\^). (19) 

Stripes are solutions satisfying zg = 0; rhombs correspond to zi = zg; squares cor- 
respond to zi — 2:^/2- Comparison of ( |l9|) with ( |l5|) restricted to the subspace 
z = (zi,0, ze,0) confirms the relationship 03 — h{9^)^ etc. 

To compute 7 and 62 in equation (16) we consider a three- mode solution 

w^,^ = ^le*"^!"- + Z2e'''="" + zse*^" + c.c, (20) 

whereki =gc(l,0),k2 = ^^(-1/2, \/3/2), and kg = (7c(-l/2, -\/3/2). The amplitudes 
evolve according to equations of the form 

ii =Ai^i+7z*z* +zi(6i|zi|2-h 62(1^21' + k3p)) (21) 

with similar equations for Z2 and 23. 

To compute the bifurcation equations (^^ and ( |2l| ) from equation (|l|) we first 
Taylor-expand f{u, v) and g{u, v) through cubic order in u and v: 

f{u,v) = au + bv + F2{u,v) + F3{u,v) + ■ ■ ■ 
g{u,v) = cu + dv + G2{u,v) + G3{u,v) + ■ ■ ■ . 



Here F2{u, w), G2(u, w), F^^u, v), and G3(m, u) denote the quadratic and cubic nonlinear 
terms: 

/ F2iu,v) \ ^ ( fuuU^/2 + fuvUV + /™ 1^2/2 \ 

/' F3{u,v)\ ^ f fuuuU^/G + fuuvU'^v/2 + fuvvUv'^/2 + fyyyV^/6\ 

\G3{u,v) J ~ \guuuU^/Q + guuvu'^v/2 + guvvUV^/2 + gyy^v^/6 J ' 

where the subscripts on / and g denote partial derivatives evaluated at u = w = 0. 

We use perturbation theory to compute the small amplitude, slow evolution on the 
center manifold. Let 

ti = et, t2 = e^t, 

^ wi{r,h,t2) + e^W2 + e^Ws + 0{e^), (22) 

where Wi takes the form of cither ( [l^ ) or (|20| ) for the respective rhombic or hexagonal 
calculations, and W2 and W3 represent higher order terms to be determined. The 
vector (ui, vi)'^ is a null-vector of L{qc), at the bifurcation point; see equation (||). We 
further assume that the linear parameters a, 6, c, d, K are within O(e^) of the critical 
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bifurcation values, e.g. a = ac + e^Xa, etc. where ac,bc, Cc, dc, Kc satisfy (||) and 
determine qc via substitution into (^. With these assumptions, equation ([^) may be 
written as 

,j 3 /u\^a fn\^, /A -.2, f u\ ^ / F,(u.v)\ ^ / 



phis higher order terms, where 



■^0=1 „ Y72 



and 



Xa 

Xc Xd + A/fV^ 



3.1 Rhombic Lattice Computation 

For the rhombic case, the equation for W2 is given at O(e^) as 

"0 tI-^i^ + ^oW, = ( ^f'^'^'DwX, (24) 
vi ) dh \G2{ui,vi) ) 

where Wrh is given by (|8|) . Through the term w^^^ the quadratic nonhnearities generate 
four terms with distinct wave numbers: 

= z2e*^''i-'' + z^e*2k«-r + c.c., V2w2 = -4g>2, 
W3 = ziZee*(''i+'"')-'' + c.c., V^wg - -29,^(1 + cos 61)^3, 
u;4 = ziZge^C'i"''''^-'" + c.c, V2u;4 = -2g2(l _ cos 6')ti;4, 



There are no terms on the right hand side of (|^) with wave number (jc, so the solvabihty 
condition is 

— Wrh = 0. 

oil 

Since there is no dependence on the ii-time scale it will be omitted from the remaining 
analysis. Solving equation (^) determines 



W2 



E :;; (25) 



where u„ and w„ (n = 2, . . . , 5) are constant and satisfy 

""^ w„ === ^ ^2(ui,wi) 
/ " VG2(mi,wi) 



r I F2(ui,Ui) , , . 

^0 = ) N c„w„, (26) 
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where C2 — 1 and cs — C4 — C5 ^ 2. Explicit expressions for it„ and Vn are given in the 
Appendix. 

The equation at 0{e^) is 

Wr'^^rh + L0W3 = L2 f Wrh + I fir ft"] „ N I U2Wrh 



Vl J dt2 V ^^1 



' aF2(«i,^^i) 

0G2("i,^'i) 
V dvi 

F3{ui,Vi 




(27) 



where {U2, V2) are given by (|25|). There are now terms on the right hand side of the 
equation with spatial dependence e**'^ '' and e*'''' ""; the operator Lq is not invertible for 
terms with this Fourier dependence. To ensure that these terms will lie in the range of 
Lq we apply the Fredholm alternative theorem to obtain the solvability condition. In 
this computation we use the fact that the nuUspace of the adjoint linear operator Lq 
is spanned by 

(^^,^)e*'^l•^ (u, ^)e*'^''■^ 

where 

{u, i) = (c, -fle + g?) • (28) 

The inner product of equation ( p7| ) with these left null- vectors yields the solvability 
condition, which takes the expected form of the bifurcation equations (pj[). After 
rescaling time to absorb the factor uiu + viv, which is always positive, we obtain the 
rhombic bifurcation equation coefScients 



and 



where 



ai = (u2 + U5)rii + {v2 + v^)rj2 + 3/3 (29) 



h[e) = (w3 + 1*4 + uz)Vi + [v^ +Vi + V5)r,2 + 6/3, (30) 



.dF2{ui,vi) ^dG2{ui,vi) 

m = U H 'rV , 

oui aui 
^dF2{ui,vi) ^dG2{ui,vi) 

(3 = uF3{ui,vi) + vG3{ui,vi). 
The expressions for m„ and w„, n — 2, ... ,5, are given in the Appendix. 
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3.2 Hexagonal Lattice Computation 

The hexagonal calculation proceeds in a fashion similar to the rhombic case, but re- 
quires some additional work due to the presence of resonant terms with wave number 
Qc which are generated by the quadratic nonlinearities. 

To simplify the presentation we let zi = Z2 = — in (|2^). Since the form of 
the amplitude equation is known from equation (p^), and 61 (= ai) is known from the 
rhombic lattice calculation, the remaining coefficients 7 and 62 may be extracted from 
the quadratic and cubic coefficient expressions, respectively. 

At O(e^) the equation to be solved is 

vi ) dti \G2[ui,vi) y ns'^ 

where Whcx is given by (po|). The term w^^^^ generates four terms with distinct wave 
numbers: 



W2 



W5 = \zh\'^, V^tus = 0, 



We 



= zr(e*''^''' + e'''^-"- + e'''^-'') + c.c., V^wy - -g>7- 

From W7 we see that there are now terms on the right hand side of ( ^2|) with wave 
number Qc, for which Lq is not invertible. Unlike in the rhombic case the ti time 
scale is necessary here, reflecting that the small-amplitude dynamics are dominated by 
the quadratic nonlinearities. As in the rhombic calculation the Fredholm alternative 
theorem is applied, giving the solvability condition 

d 2 

{uiu + viv) — Zh = jzl , j = 2{uF2{ui,vi)+vG2{ui,vi)). (33) 

Oti 

The factor {uiu+viv) is positive (from conditions (^) and (||)) and will later be absorbed 
by a rescaling of time. Note that if 7 is 0{e) then the ii-time scale is unnecessary. This 
is the situation of interest here, since we know from Table ^ that when 7 is 0(1) there 
are no stable small-amplitude steady states. However, in the following computation we 
must retain both time scales to compute the general form of the cubic coefficient 62. 
In our subsequent analysis we will focus on the degenerate problem 7 = 0. 
Continuing as in the rhombic case, the solution of ( p2| ) is 

^ ^ n=2,5,6,7 ^ 

where U2, V2, 1*5, and 115 are known from the rhombic lattice computation, and uq and 
are computed in a similar fashion. Equations (|3^) and (^) lead to the following 
solvable equation for (uj, wy): 

T f 'U7\ 7 f ui\ f F2(ui,Vl)\ 

-t^O ( . ]WY = — z ( , W7 + 2 „ { W7. 



V7 J UiU + ViV\ViJ \Ct2(Mi,Wi) 
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We choose the particular sohition of this singular problem to be 

Zi7 = 0, 1,7 = -If ;Z!^ + 2F2(ui,i;i) 

be \ uiu + viv 

At 0{e^) the equation is 

^ ^ U'hcx + 4- W2 + L0W3 = i2 f "0 ^«hc 



vi J dt2 dti \vi 

dF2{ui,vi) 

dui 
dG2(ui,vi) 
dui 

8^2(111 -vi) 

dvi 
dG2{ui,vi) 



\ dvi 
' F3{ui,Vi) 




3 

cx' 



(34) 



Again the Fredholm alternative theorem must be applied, giving a solvability condition 
of the form 

d 

{uiu + viv)—Zh = Xzh + {bi + 2h2)\zh\^Zh. 
at2 

By letting z(t) — ezh{ti, 12), /i — e^A, and rescaling time to absorb the factor uiu + viv 
on the left hand side, we obtain the reconstituted bifurcation equation 

i = + 7z*^ + {bi + 2b2)\z\^z. 

The coefficients /i, 7, 61, and &2 ai'e given in the Appendix. 



4 Parameter Collapse 

Naively one might expect that the cubic coefficients a„ and 6„ in the bifurcation equa- 
tions may be adjusted to take on any desired value, by varying the many nonlinear 
coefficients in the Taylor expansions of f{u,v) and g{u,v). However, the expres- 
sions for a„ and 6„, given in Appendix ^, clearly show that all eight cubic coeffi- 
cients contained within F3{ui,vi) and G3{ui,vi) collapse into the single expression 
P = uF3{ui, vi) + vGsiui, vi). Similarly, the six quadratic coefficients in ^2(14, v) and 
G2(m, v) are contained in the coefficient expressions solely in the four terms i^2(wi, vi), 
G2(wi,wi), rji and 772, given by (^). Moreover, these last four quantities satisfy the 
relation 

uiVi + Vim = 2{tF2(ui, Wl) + 2wG2(ui, Wl). (35) 

Therefore the six quadratic coefficients collapse into three independent parameters. In 
summary, the fourteen nonlinear coefficients in the Taylor expansions of f{u,v) and 
g{u, v) enter our analysis in just four independent combinations. 
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4.1 Hexagonal degeneracy 

The parameter collapse is even more dramatic at the degeneracy 7 = in the hexagonal 
bifurcation problem An additional striking phenomena occurs at this point, which 
results in all ^-dependence of the coefficients dropping out, thereby severely restricting 
the types of patterns which may bifurcate stably. Below we show that the bifurcation 
equations depend on just two quantities determined from the reaction kinetics of 
Thus, we may simply and completely characterize the possible bifurcation scenarios 
for all reaction-diffusion systems of the form (|l]) restricted to the hexagonal bifurcation 
problem at the degeneracy. 
The degeneracy condition is 

J = uF2{ui,vi) +vG2{ui,vi) ^ 0, (36) 

which implies 

-U1771 -I- vir]2 — (37) 

by substitution into (|35|). These two equations may be solved for F2 and rji in terms 
of G2 and i]2, and therefore these four parameters collapse into the single parameter 
combination 'q2G2{ui,vi) in the coefficient expressions for a„ and 6„. 

The bifurcation coefficients hi = oi, h{6), and &2 contain the terms it2, ^2, . . . , ue, ^^6: 
as given in the Appendix (see also equations (|29|) and (^). These terms are determined 



by solving equations of the form (c/. (26)) 



Here io('?c) is the singular matrix (^) from the linear calculation, c„ is a constant, and 
fn{d) is nonzero; for example, fz{0) = 1-1- 2cos0. Multiplying the above equation by 
the left null- vector {u,v) of Lo{qc) gives 

ilUn + KvVn — 

at the degeneracy. Using this relation and relation (p^), we see that 



/ Kviv 

UnVl + = — + 1 

\ UlU 

From the definitions of the null- vectors (ui,wi) and {u,v) it immediately follows that 

+ w„?72 =0, n = 2, . . . , 6. 

These expressions may then be substituted into the bifurcation coefficient expressions 
from the Appendix, giving 

bi = {u2 + uz)rii + {v2 + vz)ri2 + 3/? = 3/3, 

h{9) = iu^ + M4 + u<s)rii + (w3 + W4 + v<s)ri2 + 6/3 = 6/3, 
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15p+2a=0 



Hexagons 



3P+a=0 



Stripes 



Figure 2: Stability of patterns at the hexagonal degeneracy 7 = 0. The details of 
f{u,v) and g{u,v) in (|^) enter the stability analysis via two effective parameters: a, 
which depends on the quadratic nonlinearities, and /3, which depends on the cubic non- 
linearities. Rhombs and super hexagons are always unstable, as shown in Section 4.1, 
All patterns bifurcate unstably for /3 > and/or 15/3 + 2a > 0. For the region labeled 
"hexagons" we have determined only that hexagons are neutrally stable, since their 
stability for 7 = depends on terms higher than cubic in (16). 
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b2 = (u5 + U7)t]i + {v5 + ve + vr)ri2 — 2'y{vve) + 6/3 = ^7772 + 6/3, 

for 7 = 0. In short, the quadratic nonlinearities from F2{u,v) and G2{u,v) completely 
disappear from the bifurcation coefficients bi and h{0), and appear in 62 solely through 
the term 

a = vjri2. (38) 

This further implies that all ^-dependence disappears from the bifurcation equations, 
and that rhombs and super hexagons are always unstable at the bifurcation point, for 
the degenerate problem. The signs of the eigenvalues for rhombs and super hexagons 
are given in Table ^. The first two eigenvalues for rhombs depend on the signs of 61 + 64 
and hi — 64, which cannot both be negative. Hence rhombs are always unstable. The 
signs of the first three super hexagon eigenvalues are determined by the signs of 

ei = hi + 262 + ^4 + 65 + fee = 33/3 + 2a, 

62 = 61 + 262 - &4 - 65 - fee = -3/3 + 2q, 

63 = fei - fe2 = -3/3 - a. 

All three cannot be made negative, since ei = — (3e2 + 863). Stripes and hexagons 
may bifurcate stably; Figure]^ gives their stability assignments in the (a, /3)-parameter 
plane. Note that, for the cubic truncation of the bifurcation problem, hexagons are only 
known to be neutrally stable, since the eigenvalue — 7X + 0{x^) in Table g depends on 
higher-order terms, which are not computed here. For the corresponding square lattice 
computation, stripes are stable for /3 < 0; all other solutions are unstable for 7 = 0. 
This picture can be made even simpler for systems with F2{u,v) oc 6*2(1*, w), such as 
the CIMA model considered in the next section. In these systems, 7 = typically 
imphes i^2(wi,wi) = 62(1*1, wi) = 0, leading to a = in Figure |[ Stripes are then 
stable provided /3 < 0. 

We next unfold the degenerate problem by considering < I7I ^ 1 in (po|), keeping 
61 = 3/3, 62 = a + 6/3, and &4 = 65 = 65 = 6/3. There is now a small range of values 
/X e {pi, ^2) for which super hexagons (or triangles) are stable, provided that a < 0, 
/3 G {a/3, —a/21). For example. Figure || shows the unfolded bifurcation diagram for 
/3 = — 1/3, a = — 2 (0 < 7 ^ 1). It is straightforward to construct a system of the 
form (|l|) with these values. Let 

K ^A, 

f{u, v) = {2 + Xa)u - l.8v + 0.2378^2 - O.OMIm^ (39) 

g{u, v) = 4.05u - 2.8v - O.OW + uv. 

This system gives 7 « 0.03, a ~ —2, (3 w —1/3, and according to the unfolding 
calculation gives stable super hexagons (or triangles) in the range 3.07 x 10^''' < Xa < 
1.38 X lO"''. A value of = 7.6 x 10"^ corresponds to ^ = 5.54 x lO""* in Figure ||, at 
which point rolls, hexagons, and super hexagons are all stable according to our analysis. 
Figure shows a numerical integration of this system with Aq = 7.6 x 10^^, using a 
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0.0015 



Figure 3: Hexagonal lattice bifurcation diagram for the reaction-diffusion system with 
7 = 0.03, a = —2, and (3 — —1/3, as described in Section 4.1. Secondary branches and 
unstable primary branches are not plotted. Clearly there is a range of fi for which this 
system is tri-stable. 
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Figure 4: The steady state concentration field u(r) in the form of stripes, hexagons, 
and super hexagons, as computed numerically from the reaction-diffusion system (|^) 
with reaction-kinetics given by (|3S|). The sole difference between these numerical runs 
was in the initial condition. Compare with the bifurcation diagram given in Figure 
Numerical simulations were performed with a pseudo-spectral Crank-Nicholson scheme, 
using second-order Adams-Bashforth on the nonlinear terms. The rectangular domain 
has aspect ratio a/S, with a 128x128 grid. 
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pseudo-spectral Crank-Nicholson scheme on a 128 x 128 rectangular grid of aspect ratio 
All three pictures are at the same parameter values, and differ solely in the initial 
condition. The amplitude of the fundamental mode in each case is within two percent 
of the steady-state values predicted by the amplitude equations. 



5 Application to Model Equations 

5.1 Lengyel-Epstein model 

The results of the previous sections will now be applied to a reaction-diffusion system 
proposed by Lengyel and Epstein |^ as a reduced model for the CIMA reaction: 

At = a~A~^ + W^A 
Bt = 5b(A-^ + IV^B). 

Note that the nonlinear terms differ solely by a multiplicative constant, which will 
further simplify the hexagonal lattice computations. 

Our analysis assumes a steady-state solution of u = w = 0. Thus we let 

A = xa + {l+xl)u, B^{l + xl){l + v) 

where 

Xo EE - > 0. 
5 

After expanding the nonlinearities to cubic order the equation may be written as 
Ut — y'^u + au + bv + F2[u,v) + Fz[u,v) 

Vt = 5b[KV'^v + CM + + G2 (u, w) + G3 (u, w)] , ' 

after rescaling time and space to absorb a factor of 1 + Xq . The coefficients in (^l|) are 
related to those in (^0|) by 

a — Sccq — 5, b — —4x0, c = 2a;Q, d — — xq, K — - 

b 

F2{u, v) = 4xo(3 - Xo)u^ + 4(xo - l)uv, 
Fsiu, v) = 4:{x^ - 6x1 + 1)"^ + 4a;o(3 - xl)u^v, 

G2{u,v) = ^F2{u,v), G3{u,v) = ^F3(u,v). 

The overall factor of 6b in the Ui-equation will later be removed from the bifurcation 
problem by rescaling time. That is, although the overall factor of 6b enters the linear 
calculations - in particular the necessary conditions for Turing pattern formation - 
it simply scales out of the final bifurcation problem. In addition to the two remaining 
equation parameters, a and K, the lattice angle will enter our bifurcation analysis. 
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# Rolls 

# Rhombs (Rh) 

Squares (Sq) 

Rh+Sq 

Super/Anti- 
Squares 

# Unstable 



■n/2 



Figure 5: Square lattice bifurcation results (for various lattice angles 9) as a function of 
the parameter a in the Lengyel-Epstein CIMA model ( |40|) . The diagram was generated 
by substituting the coefficient expressions from Appendix ^ into Table ^ and evaluating 
at different points in the a — 9 plane. The region near lattice angle 9 = ^/3 has been 
removed, since this angle gives a hexagonal interaction,. Note also that at a « 9.4 
rolls are stable for all 9. This value of a yields the degeneracy 7 = in the hexagonal 
bifurcation problem ([l6|). 
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Choosing K to fix the system at the Turing bifurcation point, the resulting bifurcation 
scenario is then determined by a and 9, where, for example, 9 is related to the choice 
of square lattice by (0). 

We first consider the relative stability of patterns which lie on a square lattice. 
Figure]^ is a plot in the (0, d) plane, generated by calculating the bifurcation coefhcients 
oi, 12 = h{'!T/2), as = h{9), and 04 = /i(7r/2 — 9), as given in the Appendix. These 
coefficients are then substituted into the eigenvalue expressions in Table 0, generating 
the stability assignments in the figure. The region surrounding = 7r/3 has been 
removed since h{9) diverges as 6* — > 7r/3; the hexagonal lattice analysis is required at 
this point. Also, in the vicinity of 9 — tt/3 the domain of validity of the bifurcation 
results is very small; certain slaved modes are only weakly damped, leading to a small- 
divisor problem in the computation of h{9) (see Figure ^for an example). Finally, we 
emphasize that anti- and super square patterns only exist for a discrete set of 9 values, 
as discussed in Section B. 




Figure 6: Critical modes (2, 1) and (—2, 1) on a (2, l)-square lattice interact nonlinearly 
to generate a mode (0,2) that lies very near to the critical circle. We expect that this 
proximity will greatly reduce the domain of validity of the bifurcation calculation. 



A number of general conclusions may be drawn from Figure At all values of d, 
simple squares are unstable for some range of 6'-values. We therefore conclude that 
they are unstable, near onset, in the unbounded domain. Rolls are similarly unstable 
for all a with the exception of d « 9.4, where rolls are stable to perturbations in any 
direction 9. Below we discuss the significance of this particular d-value. 

Figure J^shows numerical simulations at two points in Figure |^. A square box size 
of L = ^/l3(2n/qc) gives critical lattice vectors of (±3, 2)(27r/i) and (±2, 3)(27r/i) 
(c/. Figure A (3, 2)-lattice generates mode angles of 6*3 — cos^^ 12/13 ~ 7r/8 and 
94 = 7r/2 — 03 « 37r/8. Figure ^ suggests that for 9 — Sir/ 8 rolls become unstable to 
rhombs at d « 9.8. The rolls in Figure]^ were generated from random initial conditions 
at d = 9.74; the rhombs were generated from random initial conditions at d = 9.8. For 
these simulations (5 = 1, 6 = 5, and the numerics were performed quite close to the 
bifurcation point, using a value of K = + 2.9 x 10"**. These numerical results are 
in fact sensitive to the size and shape of the computational domain; the choice of the 
computational domain damps modes that might otherwise destabilize the rhombs. 

In Figure |^, rolls appear to become stable for all 9 at the point a « 9.4. This is the 
point at which the hexagonal degeneracy condition 7 = is satisfied; since the nonlinear 
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Figure 7: Rolls and Rhombs in Lengyel-Epstein CIMA model (^) for 5 = 1, 6 = 5, and 
c chosen close to the Turing bifurcation point, as discussed in the text. Plot is of the 
numerically computed amplitude of the activator variable u, starting with a random 
initial condition. Rolls were generated at a = 9.74, rhombs at a = 9.8, in agreement 
with Figure |^. The computational domain corresponds to a (3,2)-square lattice; the 
grid size is 32 x 32. 



X 



Hex 



Stripes 



Figure 8: Unfolding of the degeneracy (a — 9.39 . . .) for the Lengyel-Epstcin CIMA 
model. All other primary solution branches on the hexagonal lattices are unstable. 
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Figure 9: Examples of super squares and anti-squares for (m,n) — (5,2), i.e. a box 
size of V29(27r/gc)- Patterns consist of four equal-amplitude critical Fourier modes, as 
described in section |[ 

terms differ by a constant, the only way to satisfy 7 = uF2{ui,vi) + uG2(ui,ui) = 
in this model is to take F2{ui,vi) = G2(mi,ui) = 0, which occurs at 



9.39. 



At this point, hexagons are unstable to rolls (c/. Figure ^ for a — 0, P < 0). 

The results summarized by Table |^ may be used to determine certain aspects of 



the hexagonal bifurcation problems near the degeneracy 7 = 0. In section 4.1 



it was 



pointed out that the unfolding scenario is fairly simple; in this model, the unfolding 
computation is even s imp ler. Since i<2(wi, vi) =0, the parameter a is zero. And since 
we know from section O that there is no ^-dependence, the unfolding is described by 
the single bifurcation diagram given in Figure H. 



5.2 An Example of Super Squares 

Section ^ described two types of superlattice patterns for square lattices: super squares 
and anti-squares. Figure ^ gives an example of each of these patterns for (m, n) — (5,2), 
corresponding to 63 = cos~^ (20/29). Each pattern has a 7r/2 rotational symmetry. 
Super squares are invariant under reflections of the square domain; anti-squares have 
symmetries involving both translations and reflections (see Q for details). 

In Section |^ a system was constructed which supported stable super hexagons; by 
using the Appendix coefficient expressions, coupled with the stability information in 
Table we may similarly construct a system which supports either super or anti- 



26 




Figure 10: Numerical integration of system ( p^ ) described in Section |5.2| , with a box of 
size •\/29(27r/(7c), and on a 64 x 64 grid. The initial (left) and final (right) distributions 
of the field u{r) are shown. Compare the final state with the super squares of Figure ||. 

squares. (The relative stability of super and anti-squares is determined by terms of 
higher than cubic order in (p^), which we have not computed.) The reaction-diffusion 
system (^ with 

f{u,v) — au — 1.8w — O.lSuw + ti^, 

g{u,v) = 6u-2.8v + 2.8v^ + 0.2u^, (42) 
K ^ 9, 

undergoes a Turing bifurcation at Oc = 1.8797 . . . with critical wave number Qc — 
0.88563 . . .. Figure |l^ shows a numerical integration of this system at a = 1.88, just 
above onset, showing both the initial condition and final steady state. 

6 Conclusions 

From the general two-component reaction-diffusion system (|l|) , we have derived analytic 
expressions for the coefficients of the leading nonlinear terms in the bifurcation equa- 
tions for rhombic, square and hexagonal lattices. These coefficients allow us to calculate 
the relative stability of patterns which are periodic on either square lattices (stripes, 
squares, rhombs, and super squares or anti-squares) or hexagonal lattices (stripes, sim- 
ple hexagons, rhombs, super hexagons or super triangles) at the onset of the Turing 
instability. In particular, the dependence of stability upon system parameters may 
be calculated explicitly for specific reaction-diffusion models, as we demonstrated for 
several model equations, including the Lengyel-Epstein CIMA model. 
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One of the surprising results of our general analysis is that details of the reaction- 
kinetics enter the computation of the bifurcation cocfhcients in a very limited fashion. 
In the case of the degenerate hexagonal bifurcation problem we show that the coeffi- 
cients of the cubic terms in the bifurcation problem depend on the reaction-kinetics 
through just two effective parameters, that are simple to compute. Moreover, we find 
that at this degeneracy the angle dependence, which is expected for the rhombic lattice 
problems, drops out completely, and that rhombs and super hexagons never bifurcate 
stably at the Turing bifurcation point of two-component reaction-diffusion models. 
However, in the unfolding of this degenerate bifurcation problem we find that super 
hexagons may co-exist stably with simple hexagons and stripes for a limited range of 
the bifurcation parameter. This tri-stability property is demonstrated by numerical 
simulation of a reaction-diffusion system in the vicinity of the degenerate bifurcation 
point. 

We plan to carry out more extensive numerical investigations that test some of the 
predictions of our bifurcation analysis. This will enable us to investigate the domain 
of validity of our bifurcation analysis, which is restricted to small-amplitude spatially- 
periodic Turing patterns. We also intend to determine the nature of the curious behav- 
ior at the hexagonal degeneracy, and to determine to what extent parameter collapse 
is exhibited by systems with three or more components, e.g. for Turing patterns in 
activator-inhibitor- immobilizer models piJl . 
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A Appendix 



Here we summarize our calculations of the coefficients in the bifurcation problems 
associated with rhombs and hexagons, for systems of the form 



Ut ~ V u + au + bv + F2{u, v) + F3{u, v) + ■ 

I _ _ _ Y72 _ a I :, 

-'yy ■ 



vt = KV'^v CU + dv + G2{u,v) + G^{u,v) ^ , V'^ = dxx + dy 



The calculation assumes that the linear coefficients are near the Turing bifurcation 
point, i.e., 

a = ac + Xa, b^bc + Xb, c = Cc + Xc, d ^ dc + Xd, K^Kc + Xk, 

{K^ac - dcf + AKcbcCc = 0, 

and that the system is in the Turing regime (conditions (^). The following values are 
common to all coefficient expressions: the right and left null vectors from the linear 
problem, 

ui \ ^ f -be \ f u\ _ f c, 
vi J ~ \ac- 

the critical wave number Qc, where 



2 KqQ,q 



2Kc 

the quadratic and cubic nonlinearities from the Taylor expansion of the reaction terms, 
G2{u,v) J ~ \guuu'^/2 + guvUv + gyyv'^/2^ 

F3{u,v)\ ^ f fuuuU^/G + fuuvu'^v/2 + fuvvUv'^/'2 + fvvvV^/6' 

Gz{u,v) ) ~ \guuuU^/6 + guuvu'^v/2 + g^yyUV^/2 + gyyyV^/6 J ' 
and three effective parameters involving the quadratic and cubic nonlinearities, 

^dF2{ui,vi) ^dG2{ui,vi) 

'^1 = " 5 5 

oui oui 

^dF2{ui,vi) ^dG2{ui,vi) 

m = u — h V — 

OVi OVi 

13 = uF3{ui,vi) +vG'i{ui,vi). 
The value of the linear term in the amplitude equations below is given by 



Finally, an overall factor of Miu + uiw, which is always positive, has been removed from 
the equations by rescaling time. 
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A.l Rhombic lattice 

The amplitude equation for the rhombic lattice problem is 

ii = fizi + ai\zi\'^zi + h{6)\z0\'^zi, 
zg = jizg + ai\ze\'^ze + h{9)\zi\'^Z0. 

The coefficients are given by 

ai = (w2 + U5)t]i + {v2 + V5)r]2 + 3(3, 

h{6) = {U3 + ^4 + U5)r]i + {V3 +V4 + V5)t]2 + 6(3, 



where 



with 



M2 \ 1 / 4Keql -de be \ ( F2{ui,Vi) 

V2 J 9Keqi\ Ce 4q'^ - tte J \G2{ui, Vi) 

U3\ Ke{ae + 1qlcos6) be \( F2{uuVi) 

vaj e+ V Cc ^^+2qlcose J \G2{ui,vi) 

Ui\ ^ 2_ f Ke{ae-2q'^cose) be \ f F2iui,Vi) 

ViJ V ^^-2q^^cos9j \G2iui,Vi) J ' 

Uo \ ^ 2 /-de \ f F2{ui,vi)\ 

V5 J Keqt\Ce -ae J \G2{ui,Vi) J ' 

C± = ^cgc(l±2cos^)2. 



It may be shown that 

U3 + U4 + U5 = — — —-r[F2{ui,vi)(2Keac-dc{l + 16cos^e)) 

Keqe{i — 4cos^ o*) V 

+G2(wi,Ui)6e(3 + 16cos''6») 

1^3 + ^4 + 1^5 = A,. ^. ■^2(wi,Ul)Cc(3 + 16C0S^6I) 

+G2(«i,vi) ('2-^-ae(l + 16cos4 0) ) ). 



A. 2 Hexagonal lattice 

The amplitude equation for the hexagonal lattice problem is 

ii = iizi + -iz^zl + bi_\zifzi + 62(1^2!^ + \z::i\^)zi, 
with Z2 and is obtained by cyclically permuting zi, Z2,Z3. The coefficients are given 

by 

7 = 2{uF2{ui,vi) + vG2{ui,vi)), 
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bi = ai, 

b2 = {U5 + We)??! + (^^5 + ^6 + V7)V2 —, Z + 6/3, 

UlU + ViV 

where ui, vi, U5, and V5 are given in the rhombic case, and 

/ \ 1 ( 2,K^ql -de be \( ^2(^1, \ 

- -2f2K,^i)V 

be \uiu + viv J 
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